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Abstract 



We study the density of the supremum of a strictly stable Levy process. As was proved recently 
in [23] , for almost all irrational values of the stability parameter a this density can be represented 
by an absolutely convergent series. We show that this result is not valid for all irrational values of a: 
we construct a dense uncountable set of irrational numbers a for which the series does not converge 
absolutely. Our second goal is to investigate in more detail the important case when a is rational. 
We derive an explicit formula for the Mellin transform of the supremum, which is given in terms of 
Gamma function and dilogarithm. In order to illustrate the usefulness of these results we perform 
several numerical experiments and discuss their implications. Finally, we state some interesting 
connections that this problem has to other areas of Mathematics and Mathematical Physics, such 
as q-series, Diophantine approximations and quantum dilogarithms, and we also suggest several 
open problems. 
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1 Introduction 



Let X be a strictly stable Levy process, which is started at zero and is described by the stability 
parameter a G (0,1) U (1,2) and the skewness parameter (3 G [—1,1]. We assume that \X\ is not a 
subordinator, which implies that |/3| ^ 1 if a G (0, 1). It is well-known that the characteristic exponent 
^(z) := — ln(E[exp(LzXi)]) is given by 

= c\z\ a (l - i/3 tan (^j sign(z)) , z G R, (1) 

where c is a positive constant (see Theorem 14.15 in [35] or Chapter 8 in [4]). The characteristic exponent 
y&(z) satisfies ^/(a^z) = a^f(z) for all a > 0, which implies the important scaling property: the processes 
{X at : t > 0} and {a^X t : t > 0} have the same distribution. In particular, without loss of generality 
we can assume that the parameter c in (1) is normalized so that 

c 2 (l + /3 2 tan(^) 2 )=l. (2) 

As we will see later, it is more natural to parametrize stable processes by parameters (a,p) instead 
of (a, (3), where the positivity parameter p := F(X 1 > 0) can be expressed in terms of (a, (3) as follows 

see Section 2.6 in [40]. Using Theorem 14.15 in [35] one can check that with our normalization (2) the 
Levy measure of X is given by 



n(dar) 



\x\ 



-l-a 



2sin(^) 



sin(7rap)l{ z>0 } + sin(7ra(l - p))l {x<0} 



dx. (4) 



It is also easy to see that under transformation (3) the set of parameters (a, (3) that we have described 
above is bijectively mapped onto a set 

A := {a G (0, 1), p G (0, 1)} U {a G (1, 2), p G [1 - a' 1 , a' 1 ]}, 

which is shown on Figure 1. We call the set A an admissible set of parameters. Note that formula (4) 
implies that when a G (1,2) and p = 1 — a^ 1 {p = a^ 1 } the process X is spectrally-positive {resp. 
spectrally-negative}. 

We define the supremum of the process X as S t := sup{X u : < u < t}. Note that due to the 

scaling property of stable processes, we have St — £«5i, thus it is sufficient to study the distribution of 
Si. Our main object of interest is the density of Si, which will be denoted by p(x; a, p) (or by p(x) if we 
do not need to stress the dependence on parameters). 

The study of the distribution of Si and its various integral transforms is an interesting problem that 
has a long history. As the following identity shows, an equivalent problem is to study the distribution of 
the first passage time := inf{t > : X t > h} for any positive h 



(r+ <t)= P(S t > h) = P (Si > hf 



Another closely related problem is the investigation of the Wiener- Hopf factorization for stable processes. 
The positive Wiener-Hopf factor is defined as 

(f)(z) = <f){z; a, p) := E [e^w] , Re(z) > 0, 
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where e(q) is independent of X and has exponential distribution with E[e(g)] — 1/q. We see that the 
Wiener- Hopf factor gives us information about the supremum evaluated at an exponential time e(g), while 
we would like to obtain the infomation about the supremum at the deterministic time. The following 
identity, which is based on the scaling property, connects these two problems 

oo 

J z w - x E [e- zS ^] dz = I»E [(Se(i)P] = r M E [e(l) _ - (Si)~ W ] = I»r (l - ^) E [(S^] . (5) 
o 

The second piece of the puzzle is provided by the integral representation for the Wiener-Hopf factor 

In (E \e-* s **]) = — /in ( q — ) Re(z) > 0. (6) 

K 1 iJ 2vr J \q + V(u)J u{u-\zY KJ KJ 

Formula (6) is in fact valid for all Levy processes satisfying a mild regularity condition (see Lemma 4.2 
in [29]). 

It seems that Darling was the first to investigate the supremum of a stable process. In the paper 
[10] which was published in 1956 he has obtained (5) and (6) for symmetric stable processes (p = 1/2). 
As an application of his results Darling has found a simple expression for the density of the supremum 
of Cauchy process ((a, p) = (1,1/2)). Then in 1969 Heyde [22] has shown that formulas (5) and (6) 
are also valid for general non-symmetric stable processes. Many interesting results were discovered by 
Bingham and were published in his 1973 paper [6]. Bingham has also used an approach based on (5) 
and (6) and has obtained an absolutely convergent series representation for the density of Si when the 
process X is spectrally negative (in this case the distribution of S\ is related to Mittag-Leffier functions). 
Note that the spectrally-negative case is rather simple, as it follows from the general theory that the 
supremum S e ( q ) for any spectrally-negative Levy process must necessarily have exponential distribution 
(see formula (8.2) in [27]). In the general case Bingham has obtained the first asymptotic term for 
P(Si < x) as x — > + . His results also provide the first indication of the fact that the analytic structure 
of the distribution of the supremum may be related to certain arithmetic properties of the parameters: 
he derives a general asymptotic result for the Wiener-Hopf factor (j)(z), and this result is valid only if a 
is irrational (see formula (9) in [6]). 

Bingham [6] states that in general Darling's integral (the integral in the right-hand side of (6)) 
can not be evaluated explicitly, and that the approach based on (5) and (6) can only succeed in the 
spectrally-negative case. That this prediction was not entirely correct was shown in 1987 by Doney [11]. 
Among other results, Doney has shown that Darling's integral can be evaluated explicitly if the process 
X belongs to the class Cj-j, which is defined by the condition 

p + k = -, (7) 
a 

where k and I are integers. These classes of processes can be considered as generalizations of the spectrally 
one-sided processes, for which p — 1/a or p — 1 — 1/a. It is easy to see that the curves p = {I /a} (here 
{x} G [0, 1) denotes the fractional part of i 6 1) are dense in the set of admissible parameters A. 
These curves in fact create a rather interesting pattern (see Figure 1), and as we will see later, the 
dark/white regions on this graph correspond to different series representations for the density p(x; ot,p). 
The results by Doney clearly illustrate the importance of the arithmetic properties of parameters. We 
see that depending on whether a is rational or not the formula (7) defines a countable/dense or a finite 
set of p for which there exists an explicit formula for the Wiener-Hopf factor <p(z). 
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0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 0.2 0.4 0.6 0.8 1 1.2 1.4 1.6 1.8 2 

(a) N = 10, a e (0, 2) and p G (0, 1) (b) N = 100, a e (0, 2) and p £ (0, 1) 




0.5 0.55 0.6 0.65 0.7 0.75 0.8 0.85 0.9 0.95 1 1.5 1.55 1.6 1.65 1.7 1.75 1.8 1.85 1.9 1.95 2 



(c) N = 100, a £ (1/2, 1) and p £ (1/4, 3/4) (d) N = 100, a £ (3/2, 2) and p £ (1/3, 2/3) 

Figure 1: The set of admissible parameters A and the curves p = {la~ 1 } for \l\ < N. 



In the 1990s and 2000s there have appeared many other important results related to extrema of stable 
processes. Bertoin [5] has computed the Laplace transform of the first exit time from a finite interval 
in the spectrally one-sided case, the result is given in terms of Mittag-Leffier functions. An interesting 
duality result relating space-time Wiener-Hopf factorization of stable processes with parameters (a, p) 
and (l/a,ap) was discovered by Fourati [17], a special case of this duality was given earlier by Doney 
(see Theorem 3 in [11]). An absolutely convergent series representation for p(x) was obtained by Bernyk, 
Dalang and Peskir [3] in the spectrally-positive case. Doney [12] has found the first asymptotic term 
of p{x) as x — > +oo in the spectrally-positive case, and Patie [32] has given a complete asymptotic 
expansion. The spectrally positive case and various connections with Mittag-Leffier functions were 
investigated by Simon [36]. In the general case the first term of the asymptotic expansion of p(x) as 
x — > + or x — > +oo was obtained by Doney and Savov [13]. Peskir [33] and Simon [37] study the first 
hitting time of spectrally-positive stable processes. Graczyk and Jakubowski [19] (and independently 
Kuznetsov [25]) have discovered a series representation for the logarithm of the Wiener-Hopf factor: 
as we have mentioned above, an asymptotic version of this result was found earlier by Bingham (see 
formula (9) in [6]). The first passage time of stable processes was studied in a recent paper by Graczyk 
and Jakubowski [20]. 

In [25] we have followed the approach based on (5) and (6), and we have obtained explicit formulas for 
the Wiener-Hopf factor <p(z) and for the Mellin transform of S\ in the general case. We have also found a 
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complete asymptotic expansion for p(x) as x — > + or x — > +00, and in the case when X G Ck,i we have 
proved that this asymptotic expansion gives an absolutely convergent series representation for p(x). The 
key ingredient which allowed us to derive all these results is the so-called double gamma function G(z; r) 
(see Section 3), which was introduced and studied by Barnes in [1], [2]. Hubalek and Kuznetsov in [23] 
have proved that the asymptotic series for p{x) obtained in [25] in fact provides an absolutely convergent 
series representation for almost all irrational values of a. In this result the arithmetic properties of a are 
crucial: The result is valid for all irrational a which are not "too close" to rational numbers. Kuznetsov 
and Pardo in [26] have shown that it is also possible to derive some of the results in [25] as well as 
several new results (such as the entrance law of the excursion measure of a stable process reflected at the 
supremum/infimum) using a different technique, which is not based on (5) and (6). In some sense this 
approach is more direct and more probabilistic in spirit, it's main ingredients are exponential functionals 
of hypergeometric Levy processes and the Lamperti's transformation for positive self-similar Markov 
processes and on . 

In this paper we have three main objectives. First of all, the proof of the convergence of the series for 
p(x) given in [23] leaves open a possibility that this result is not only true for almost all irrational a, but 
may be valid for all irrational a. As we show in Section 2, this conjecture is not true: we construct an 
uncountable dense set of irrational a for which the series does not converge absolutely for almost all p. 
Our second goal is to investigate the case when a is rational: in this case the infinite series representation 
for p(x) given in [25] and [23] is not well-defined. This case is clearly important for applications, after all 
when we perform numerical computations we can only work with rational approximations to irrational 
numbers. In Section 3 we will show that when a is rational the Mellin transform of Si can be given 
explicitly in terms of rather simple and easily computable functions (Gamma function and dilogarithm) . 
Our third goal is to compute the density of Si numerically. Given the unusual properties of various 
representations of p(x), it is a non-trivial question of whether these results can be used for numerical 
computations of p(x). As we will see in Section 4, the answer to this question is positive, and the 
density of the supremum can be computed quite easily both in the case when a is rational and when 
a is irrational. Finally, in Section 5 we discuss the connections that this problem has to other areas of 
Mathematics and Mathematical Physics, and we also suggest several open problems. 

2 Convergence properties of the series representation of p(x) 

Let us first introduce several notations and definitions. For x G R we define the "floor" function 
[x\ := max{n G Z : n < x}, the fractional part {x} := x — \_x\ and we let ||x|| := min{|x — n| : n G Z} 
denote the distance to the nearest integer. For any real a we define the set 

C(a) := {p G (0, 1) : p = {la' 1 } for some / G Z}. (8) 

Thus p G C(a) if and only if X G C^i for some integers k and /. Note that if a = m/n for some co-prime 
m and n, then C(a) is the finite set 

C(m/n) = {j/m : j = 1, 2, . . . , m — 1}. 

The sets C(m/n) are clearly visible on Figure 1, they consist of points of intersections of different curves 
p = {l/a}. If a is irrational then the set C(a) is countable and dense in (0, 1). 
The following set of real numbers was introduced in [25] and [23]: 
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Definition 1. Let C be the set of all real irrational numbers x, for which there exists a constant b > 1 
such that the inequality 



x 



P 
Q 



< 



b" 



(9) 



is satisfied for infinitely many integers p and q. 



As was proved in [23], the set C is rather small: it has Hausdorff dimension zero, which implies 
that it has Lebesgue measure zero. This set is closed under addition/multiplication by non-zero ratio- 
nal numbers, therefore it is dense. It can also be conveniently characterized using continued fraction 
representation of a real number, which is defined as 



x = [a ; a 1 , a 2 , . . . } = a + 



a 1 + 



a 2 + 



where Oq G Z and Oj 6 N for % > 1. For x ^ Q the continued fraction has infinitely many terms; 
truncating it after n steps gives us a rational number p n /q n = [ao] a\, a%, a n ], which is called the n-th 
convergent. As was proved in [23], x G C if and only if there exists a constant b > 1 such that the 
inequality a n+ \ > b Qn is satisfied for infinitely many n. 

For any irrational a we define sequences {a min } m >o, n >o an d {b m>n }m>o,n>i as follows 



VI 11 + II 



r(i 



n 



— ) T(ap + m + an] 



n" 1 sin (I (ap + j — 1)) -A- sin(7ra(p + j- 1)) 



sin 



sin(7ra;j) 



r (l — p — n — —J r(o;p + m + an) 



r(i + n+^)r( 



-m — an) 



The following Theorem was proved in [23] (the statement related to the asymptotic behavior of p(x) 
has appeared earlier in [25]). 

Theorem 1. Assume that a ^ C U Q. Then for all x > 

p(x) = x-^YsY.^^"' if « £ (M), (12) 

n>0 m>0 

p(x) = x a o- l YJ2 a ^ xm+Qn > if a G (1,2). (13) 

n>0 m>0 

T/ie series converges absolutely and uniformly on compact subsets of (0, oo). Moreover, for every a ^ Q 
i/ie series (12) { (13) } provides complete asymptotic expansion as x —> +oo {resp. x —> + }. 

The above theorem, if suitably interpreted, is also true if a G £ U Q and p G C(ck). In this case 
X G Ck,i for some integers and I, and it turns out that expressions (10), (11) for coefficients a m ^ n and 
b m ^ n can be simplified, and the resulting formulas are well-defined even for rational values of a. See 
Theorem 10 in [25] for all the details. 

If we are only interested in the asymptotic expansion of p(x), then Theorem 1 and Theorem 10 in 
[25] give us a complete picture. Ifa^QjaGQ and p G C(a) } we have explicit asymptotic expansion 
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due to Theorem 1 { resp. Theorem 10 in [25] }. At the same time, in the remaining case when a G Q 
and p ^ C(a) we know that the asymptotic expansion will contain some logarithmic terms coming from 
the multiple poles of the Mellin transform of Si (see the proof of Theorem 9 in [25]). Note that even if 
a G Q and p ^ C(a) the first terms of the series (12), (13) which are well-defined still provide the first 
asymptotic terms for p(x): this follows easily from the proof of Theorem 9 in [25]. 

The situation with the absolute convergence of the series (12), (13) is much less clear. Let us 
summarize the present state of knowledge. 

(i) If a G Q then the series (12), (13) are well-defined and converge only for a finite set of p G C(a). 
This case corresponds to the points of intersections of black curves and to "white vertical intervals" 
on Figure 1. 

(ii) If a ^ Q then the series (12), (13) converge for a dense countable set of p G C(a). 

(iii) If a ^ Q and a ^ £ then the series (12), (13) converge for all p G (0, 1). 

We see that the only remaining case for which we don't know anything about the convergence of the 
series (12), (13) is when a G £. The most aesthetically pleasing and natural conjecture would be that 
condition a ^ Cm. (iii) can be dropped: 

Conjecture: if a ^ Q then the series (12), (13) converge for all p G (0, 1). 

Our main result in this section is that this conjecture is in fact false. We will show that by slightly 
modifying the definition of the set £ we can construct an uncountable dense subset £ C £ such that for 
every a 6 £ the series (12), (13) does not converge absolutely for almost all p G (0, 1). This set £ is 
defined as follows. 

Definition 2. Let £ be the set of all real irrational numbers x, for which there exist e > and b > 1 
such that the inequality 

< b-" ln ^ 1+e (14) 

is satisfied for infinitely many integers p and q. 

As the next Proposition shows, the structure of the set £ can also be conveniently described in terms 
of the continued fraction representation of real numbers. The proof of this Proposition is omitted, as it 
is essentially identical to the proof of Proposition 1 in [23] . 

Proposition 1. 

(i) If x G £ then zx G £ and z + x G £ for all z G Q \ {0}. 
(ii) x G £ if and only if x^ 1 G £. 

(iii) Let x = [do; ai, ci2, ■ ■ ■ ]■ Then x G £ if and only if there exist e > and b > 1 such that the 
inequality 

a n+1 > b qnln(qn)1+t 

is satisfied for infinitely many n. 



x 
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We see from Proposition 1 that the set C is an uncountable set, it is closed under addition and 
multiplication by non-zero rational numbers, in particular it is a dense set. It is of course a rather small 
set: being a subset of C it also has Hausdorff dimension zero and, therefore, Lebesgue measure zero. 

The next Theorem is our main result in this section. 

Theorem 2. For every a G C there exists a set of real numbers 13(a) having zero Hausdorff dimension 
(and, therefore, zero Lebesgue measure) and such that for every p £ B(a) and every x > the series 
(12), (13) do not converge absolutely. 

Proof. Assume that a G (1, 2) U C According to Definition 2 there exist b > 1 and e > such that 

\\aq\\ < qb- qln{q)1+ \ for all q G Q, (15) 

where Q is an infinite subset of N. We define the following set of real numbers 

B(a) := {p G (0, 1) : \\atp + na\\ < n ~ lnhl{1+n) for infinitely many n G N}. (16) 

From the definition (8) of the set C(a) and property (15) it is clear that C(a>) C B(a), in particular B(a) 
is non-empty. Our first goal is to prove that the set B(a) is in fact quite small: the Hausdorff dimension 
of this set is zero. 

For v > let us define 

V v (a) := {x G M : \\x + na\\ < n~ v for infinitely many n G M}. 

As was proved by Bugeaud (see Theorem 1 in [7]), for any irrational a and any v > 1 the Hausdorff 
dimension of the set V v (a) is 1 /v. From (16) we find that for any v > 

{ax : x G B(a)} C V v (a), 

and since v can be taken arbitrary large we see that the Hausdorff dimension of B(a) must be zero. 
Next, we will show that for any p ^ B(a) than for any i>0we have 

loto,?^! > 1 for all large enough q G Q, (17) 

which of course implies that the series (13) can not converge absolutely. 
First of all, from (10) we find that for n > 1 

sm(irp) T(p + n) rr , / nmTT 1 , s 

i a ».»i = ; rZ + ai) n i - i»i n ^-y, ^ 

7r (2ra)! J n\\(xn\\ 

where the first equality is obtained by applying the reflection formula for the Gamma function, and the 
inequality in the second line is derived using the following trivial estimates: 

r( P + n)>r(i) = i, 

T(ap + an) < T(l + 2n) = (2n)!, 

| sin(7T2;)| > | sin(7T2;)| < 1 and | sin(7T2)| < 7r||^||, 
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Next, since p £ B(a) we know that p C(a), in particular \\ap\\ > 0. At the same time, since 
p ^ B(a) there exists a constant C = C(a, p) > such that for all j > 1 

||ap + ja|| > Cj- lnln(1+j) . (19) 

Using (19) and the estimate 

n-l 

ln(j) In ln(l + j) < n ln(n) In ln(l + n) 

i=i 

we find 

n n—\ 

Y[\\a{p + j - = \\ap\\Y[\\ap + aj\\ (20) 

i=i i=i 

> ||a J o||C n ^ 1 exp ^-^ln(j)lnln(l+j)j > HapHC^^-™ 1 ^ 1111 ^ 1 ^. 

Next, we take q G Q, combine (15), (18), (20) with the trivial estimate (2n)\ < (2n) 2n and find that for 
every x > 

I x gi > sm(7rp)||qp|| c<q _ 1 ^ c _ 2(i , ln(2(?) _ (;ln((f)lnln(1+(;) 1 



7T 2 \\ a Q\\ 

> ^(^P) ll Q Pll ^-1^9^29 ln(2g)-gln(g)lnln(l+ g )^ln(g) 1 + e -l_ 
7T 2 

The right-hand side of the above inequality can become arbitrarily large if q is large enough. Given that 
Q is an infinite set we conclude that that for any x > it is true that |ao, ra ^ n | > 1 for infinitely many n, 
which shows that the series (13) does not converge absolutely. 

The proof in the case a G (0, 1) is very similar, the details are left to the reader. □ 



3 Mellin transform of Si when a is rational 

The fact that the Mellin transform is a powerful tool for studying stable processes was known already 
in 1950s (see Zolotarev [39] and Darling [10]). This is not a very surprising fact, given that the scaling 
property for stable processes is multiplicative in nature. We will use the following notation for the Mellin 
transform of Si 

oo 

M{w) = M{w; a, p) := E [(S^' 1 ] = J p{x; a, p)x w ~ 1 dx. (21) 



The analytical structure of this function was completely described in [25] . In particular, for all admissible 
parameters (a, p) there exists an explicit expression 

M ( s \ = a s-i G(ap;a) ^ G(a(l - p) + 2 - s; a) ^ G(a - 1 + s; a) 

{> G(a(l - p) + 1; a) G(ap-l + s;a) G(a + l-s;a)' { ' 



9 



Here G(z;t) is the so-called double gamma function, which was studied by Barnes in [1], [2]. This 
function is defined for | arg(r)| < 7T, and all z G C by a Weierstrass product 



2 



Z 2 / Z \ z 

G(z; t) := -e a r +b ^ ]~J 11+ e "s^+2(^T^ 

r m>o «>o ^ mr + nj 

m+n>0 

In particular we see that G(z; r) is an entire function of z, which has zeros at points — mr — n, n, m > 0. 
Barnes [1] has proved that it is possible to choose the coefficients a = a(r) and b = b{r) so that G(l; r) = 1 
and G(z; r) satisfies the two functional equations 

G{z + 1-t) = t(^G{z-t), G{z + T-T) = {2ny-^T- z+1 *T{z)G{z-T). (23) 

Formula (22) can be considerably simplified when X G Ckj, in this case Ai(s) can be given in terms of 
Gamma function and trigonometric functions (see formula (6.10) in [25]). As we will see in this section, 
formula (22) can also be simplified in the case when a is rational. In order to state our results, first we 
need to present several definitions. 

Everywhere in this paper we will work with the principal branch of the logarithm, which is defined 
in the domain | arg(z)| < it by requiring that ln(l) = 0. Similarly, the power function will be defined as 
z a = exp(aln(z)) in the domain | arg(^)| < ir. The principal branch of the dilogarithm function Li2(z) is 
defined in the domain z e C \ [1, oo) by the integral representation 

Ii>«:=- r^^du, (24) 
J u 



where we integrate along any path from to z which lies in C \ [1, oo) (see [28] and [31]). By expanding 
the logarithm into Taylor series it is easy to see that in the open disk {z G C : \z\ < 1} the dilogarithm 
is given by an absolutely convergent series 

Li 2 W = E^- ( 25 ) 

k>l 

While dilogarithm is not one of the elementary functions, it is a well-understood function which satisfies 
many functional identities. As we will see later in Section 4, the dilogarithm can be easily computed to 
high precision in its entire domain of definition C \ [1, oo). 

For a,q G C such that \a\ < 1 and \q\ < 1 we define the modified q-Pochhammer symbol 

ra-l 

[a;g] n :=n(l-o9 fc )». ( 26 ) 

k=l 

and for any co-prime positive integers m and n we define 



1 , ,.A (l - e 2 ™) 1 ™ 



2-nimn " v ' ) 2™ 2™ 

6 m . 6 m 



H m , n {s) := exp ( Li 2 (e 2 ^) ) x r 2ms \ mn , ' ^ ■ (27) 

e n ; e 'i 



Note that H m n (s) is well-defined and is an analytic function in the half-plane Im(s) > 0. The following 
Theorem is our main result in this Section. 
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Theorem 3. Assume that a = m/n where m and n are co-prime natural numbers. Then for all s e C 
with Im(s) > we have 

M(s) = 4 l^ e ^ m2+n2 - 3mn)+7ri (™- p ) {s -V Ei^ H m , n {mp)H m , n {ns) 

W V™ r(l-^(l- S ))^ m >(5-l)+mp)- 

Proof. We start with the following identity 

r ( i 

G(a-l + s;a) = (2tt)^V^ - G(s; a), 

which follows from (23). The above identity and formula (22) give us 

M ( ) = ( 9 \s=i -i r(s) G(ap;q) G(q(l - p) + 2 - s; a) a) 

U 1 ' r(l-^)G(a(l-p) + l;a) G(ap - 1 + s; a) G(a + 1 - s; a) ' 1 J 

The first main ingredient in our proof is the following formula 

G(z;-) = (2np n -^n^ m+n ^ +1 TT G ( — + -) TT G (- + - V' (30) 
V n / - LJ - \ m ra / - LJ - n/ 

0<fc<m-l v 7 0<fc<m-l v 7 

0<Kn-l 0<Kn-l 

k+l>0 

where G(z) := G(z; 1) is Barnes G-function. The above formula is just a special case of a more general 
result connecting G(z; — r) with G(z; r), which was established by Barnes in [2]. In order to obtain the 
identity (30) one has to set r = W\ = w% = 1 in the formulas on pages 302 and 359 in [2]. Applying 
identity (30) to each double gamma function in (29) and simplifying the resulting expression we obtain 



m — n 1 



M(s) = (27i)fa-i ^ F m , n {mp)F m , n {ns) 

r (1 - i=£) F m>n (n(s - 1) + mp) ' [6 } 



where we have defined 

m— 1 n—1 



= n n z % ^ 

k=0 1=0 \ m n nmJ 



We change the indices k H- m — 1 — k and I t— >■ n — 1 — Z in the denominator in the right-hand side of 
(32) and obtain an equivalent expression 



m— 1 n— 1 



^,n(s) = n n G^rl^Ly ^ 

k=0 1=0 \ nrn m n) 

The second main ingredient of the proof is the following reflection formula for the Barnes G-function 

- = (27T) 2 - 1 (1 - e 2 ^) 1 "^^ (^-2^|)-^Li 2 (exp(2 7 ri,)) ) ^ > Q ^ 

The above formula follows from identity (A. 5) in [38] and formulas (4.6), (8.43) in [28]. Next, we apply 
the reflection formula (34) to each term in the product (33), we simplify the resulting expression with 
the help of identities 

n— 1 „ Ji—l 



X> = ^n(n-1), ^ 2 = ^(n-l)(2n-l), 



k=l k=l 
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and we finally obtain 



F m>n (s) = (2-kY - , 

m—1 ri—1 

X 

k=0 1=0 



s- S±= „ {s~m-n)+-P-(m 2 +n 2 +3mn) 



nn(i--p(2«(-+-+ ! 

- LJ - - LJ - V V \mn m n 

k=0 1=0 

m— 1 Ti—1 
k=0 1=0 



X 



1- 



s I fc I _Z_ 
mn m n 



(35) 



Formula (35) can be further simplfied with the help of the following two identities, which are valid for 
any co-prime positive integers m and n: 



n-l 



71-1 



k=0 



n - 



l--2 n , y^Li 2 1 z< 

k=0 



n 



-U 2 {z n ) 



\z\ < 1. 



(36) 



The first identity can be easily verified if one considers the left-hand side as a product involving n roots 
of the polynomial 1 — z n . For the second identity see (3.12) in [31]. 

Using (35), (26) and (36) we can express the function F m>n (s) as follows 



jTi /„\ fr,\s~^^ ^^(s-m-n)+- p p—(m 2 +n 2 +3mn)-^— Li 2 (exp(27ris)) ^ 

Formula (28) now follows from (27) and (31), (37). 



1 s 

M g27ris^ 1 mn 



2ty\s 


2nin 




2nis 2nim 


e m 


e Tu- 




e n ■ e n 






rn 





(37) 



□ 



Theorem 3 should be compared with Theorem 2 in [25], which gives an explicit expression for the 
Wiener- Hopf factor <fi(z) in terms of the Clausen function, defined for real values of x as (^(x) := 
Im [Li 2 (exp(ix)] (see chapter 4 in [28]). This shows that when a is rational both the Wiener-Hopf 
factor <j)(z) and the Mellin transform Ai(s) can be expressed in terms of elementary functions, Gamma 
function and dilogarithm. It would be an interesting exercise to obtain expression (28) for the Mellin 
transform directly from Theorem 2 in [25] using identity (5), or, alternatively, using the approach based 
on exponential functionals and Lamperti's transformation (see [26]). 



4 Numerical experiments 

In this section we study the problem of numerically computing the density of the supremum S\. There 
are two different approaches to this problem depending on whether a is rational or not. In the case 
when a is rational we can compute p(x) via the inverse Mellin transform of Ai(s) given by (28). If a 
is irrational, the Mellin transform Ai(s) given by the general formula (28) is rather hard to evaluate 
numerically, thus we will try to compute p(x) using the series expansions (12) or (13). However, apriori 
it is not clear whether these expansions can possibly lead to any meaningful numerical results. According 
to Theorem 1 we need to ensure that a ^ C U Q, but in any computer program a would be given to a 
finite precision. Therefore, we will be working with a rational approximation to a, and we know that for 
rational values of a the coefficients of the series (12) and (13) are not even well-defined. As we will see 
in this Section, this method based on series expansions can still be used for computing p(x), though it 
certainly has its limitations. 
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4.1 Computing p(x) by inverting the Mellin transform for rational a 

In the method based on inverting the Mellin transform we face the following two problems. First, we 
need to be able to evaluate the Mellin transform Ai(s) given by (28) efficiently and with high precision. 
Second, we have to compute the inverse Mellin transform numerically. 

Let us describe how one can compute the Mellin transform A4(s). We see that the only challenging 
part in (28) is computing the dilogarithm function Li2(z) in the unit disk D := {z G C : \z\ < 1}. Let 
us denote 

Di := D n {z G C : |1 - z\ < 1/10}. 
In the domain Di the function Li 2 (z) can be computed using the identity (3.2) in [31] 

7T 2 

Li 2 (z) = -Li 2 (l - z ) + — - ln(l - z) In (z . (38) 

o 

The term Li 2 (l — z) in the right-hand side of (38) is evaluated using the series expansion (25); since 
|1 — z\ < 1/10 in Di this series converges very fast. 

In the domain D \ Di we will use the following formula 

^w---T + -K^) + -B- 1 )"^e 2 ". -mi-* (3 9) 

v 7 n>l 

where C{2n) — 1 = 2~ 2n + 3~ 2n + . . . . Formula (39) follows by combining (4.3) in [31] and the Taylor 
series of ln((27ri + w)/(2ni — w)). Let us investigate the convergence rate of (39). First of all, one can 
check that for all n > 1 we have 

2 2n (C(2n) - 1) < 4(C(2) - 1) < 3. 

Next, for all z £ ID) \ D\ we have 1/10 < |1 — z\ < 2 and | arg(l — z)\ < 7r/2, which implies that 

M 2 = I hi(l - z)\ 2 = In |1 - z\ 2 + (arg(l - z)) 2 < ln(10) 2 + tt 2 /4. 

Therefore, the n-th term in (39) can be bounded from above as follows 

C(2n) - 1 f\w\\ 2n o n f\w\\ 2n /ln(10) 2 + vr 2 /4\" , N „ 3 

^—^ — < 3 x 2" 2n ^ < 3 x v ; — '— ps 3 x 0.04919... n < . 

2n + l \2n J \2n J \ 16tt 2 J v ' 20" 

This shows that the series (39) also converges very fast. 

Once we are able to evaluate the Mellin transform (28) numerically, it is a simple matter to find p(x) 
by computing the inverse Mellin transform. Using the fact M(s) = Ai(s) we rewrite the inverse Mellin 
transform as 



p( x ) = — / M(s)x- S ds = — Re / M(l + iu)e- mHx) du 
2vri J ttx J 



1+iM 



(40) 



Thus we need to compute the Fourier transform of + iu). Note that the function A4(l + iu) is 
smooth, and as we know from Lemma 3 in [25], it decays exponentially fast as u — > +oo, thus we can 
truncate the integral in (40) at some large number (we truncate at u = 40) and then use Filon's method 
[16] to compute this Fourier integral numerically. 
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(a) Real (black) and imaginary (blue) parts of (b) The density of Si 

M{l + m) 

Figure 2: The Mellin transform and the density of the supremum in the case a = 3/2 and p = 3/5. 

The results of the computations for a = 3/2 and p = 3/5 are presented on Figure 2. Figure 2a shows 
the real and imaginary parts of the Mellin transform A4(l + iu), while on Figure 2b we plot the density 
p(x). One interesting fact that we can observe from this picture is that the density of Si is not convex. 
This fact should be compared with the general result of Rogers [34] , which states that if a Levy process 
X has a completely monotone Levy density (such as n(dx) in (4)) then for any q > the density of S e r q ) 
is also completely monotone, in particular it is convex. 

4.2 Computing p(x) using the series expansions 

In this section we investigate the second possible approach towards computing p(x) - via the series 
expansions (12), (13). The coefficients of these series are not defined for a = 3/2, thus we perturb 
a = 3/2 by a small irrational number. We will first set the parameter values as a = 3/2 + y/2/50 and 
P = 3/5. 

The implementation of (12), (13) is rather straightforward. We truncate the convergent series (13) 
at m = n = 200 and the asymptotic series (12) at m — n — 15. All computations are performed in 
Fortran90, using quad 128-bit format, which gives us working precision of approximately 34 decimal 
digits. 

We see that the results presented on Figure 3 are rather encouraging. First of all, the graph of p(x) 
on Figure 3a is very similar to the one on Figure 2b. Second, as we see on Figure 3b, in the interval 
x G (4, 5) the convergent series agrees very well with the asymptotic series. 

In order to give more credibility to this method we have performed one additional experiment. Using 
the same approach as above, we have computed p(x; a,p) for a = 3/2 ± v^2/50 and p = 3/5, and then 
we have compared the average 

p{x) := l - (p{x- 3/2 + v^/50, 3/5) + p{x; 3/2 - V2/50, 3/5)) 

with p(x;3/2,3/5) which was evaluated using the inverse Mellin transform technique described in the 
previous section. We would expect that p{x) should be close to p(x; 3/2, 3/5), since the function p(x; a, p) 
is continuous in parameter a. This continuity property can be established using the inverse Mellin 
transform representation (40) and the fact that Ai{s) is continuous in a and decays exponentially in 
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(a) (b) 

Figure 3: Computing p(x; a, p) via convergent series (13) (blue) and asymptotic series (12) (black). The 
parameters are a = 3/2 + a/2/50 and p — 3/5. 




(a) p(x; a, p) for a = 3/2 + a/2/50 (blue), (b) p(x) - p(x; 3/2, 3/5) 

a = 3/2- a/2/50 (red) and a = 3/2 (black) 

Figure 4: Computing p(x; a, p). Here a £ {3/2,3/2± a/2/50} and p = 3/5. 



|Im(s)| (with more work one can probably prove that p(x; a, p) is in fact a smooth function of all three 
variables (x,a,p)). The results of this experiment are presented in Figure 4. We see that the error 
p(x) — p{x\ 3/2, 3/5) is very small, of the order 10 -4 . This confirms that both methods are in fact 
producing reasonably accurate results. 

We would like to note that the method based on the series expansions (12), (13) seems to have some 
serious limitations. In particular, we were not able to obtain good results in our last experiment when 
we modified a = 3/2 ± a/2/100. It seems that when a is too close to a rational number, we would have to 
include more terms in the series (12), (13) and at the same time we would have to increase the working 
precision. 

We would also like to mention that there is another possible approach to computing p(x). We 
assume that a is irrational. In this case the set C(a) defined by (8) is dense, thus we can find a good 
approximation to p of the form {l/a}. This is a problem from the theory of inhomogeneous Diophantine 
approximations, and it can be solved efficiently using Cassels algorithm [9]. This procedure would give 
us parameters (at, p), for which p is close to p and p £ C(a). The process X defined by parameters 
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(a, p) belongs to Doney class Ck,i for some integers k and I, therefore we can compute p(x;a,p) using 
the simpler series expansions given in Theorem 10 in [25]. The important fact is that the coefficients of 
these series expansions involve simple bounded products of trigonometric functions. While we did not 
implement this method, we expect that in this case the computation can be performed more efficiently 
and may not even require high precision arithmetic. 



5 Conclusion and several open problems 

In conclusion we would like to discuss several connections that the study of the distribution of extrema 
of stable processes has to other areas of Mathematics and Mathematical Physics. 

First of all, the series expansions of the density of Si are related to g-series (see an excellent book 
by Gasper and Rahman [18]). The g-series can be informally defined as hypergeometric series where the 
terms involving Gamma function and factorials are replaced by the q-Pochhammer symbol, defined as 

n— 1 

(«; q) n ■= n( 1_G ^')- 

It is clear that the finite products involving sin(-) function in the definition (10) of coefficients a m ^ n can be 
rewritten in terms of the g-Pochhammer symbol, thus the series (12), (13) can be considered as g-series. 

Convergence properties of g-series when |g| = 1 is an interesting question that has attracted the 
interest of many researchers. The first results in this area were obtained by Hardy and Littlewood [21], 
where they have investigated the series 



„>« («;«)» 

Driver et. al. [14] have investigated the convergence properties of the series 

^2(a;q) n z n , 

n>0 

and Lubinsky [30] and later Buslaev [8] have studied the Rogers- Ramanuj an function 

S (?;?)» 

This function was central in Lubinsky's disproof of the Baker-Gammel- Wills conjecture (see [30]). 

The convergence properties of g-series on the unit circle |g| = 1 are naturally related to Diophantine 
approximations and properties of continued fraction expansion of arg(g). As we have seen in the proof 
of Theorem 2, in our case the situation is even more delicate and the convergence depends essentially on 
the properties of both parameters a and p. The essential tool in the proof of Theorem 2 was the result 
bu Bugeaud [7] from the theory of inhomogeneous Diophantine approximations. 

There also exists an intriguing connection between the Wiener-Hopf factorization of stable processes 
and certain special function called quantum dilogarithm, which has applications in Quantum Topology 
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and Cluster Algebra Theory. For b e C with Re (b) > the Faddeev's quantum dilogarithm is defined in 
the strip |Im(z)| < ^Re(6 + b~ l ) by 



where the singularity at x = is circled from above, see [15], [24] and the references therein. The 
function z i— > can be analytically continued to a meromorphic function on the entire complex 

plane. Combining formulas (4.6) in [24] and (4.2) in [25] we obtain a particularly simple expression for 
the positive Wiener-Hopf factor in terms of Faddeev's quantum dilogarithm 



We would like to conclude this paper by stating three open problems. 

Problem 1. Doney's Ck,i classes are well-understood from the analytical point of view. For example, 
if we consider the Wiener-Hopf factor 0(exp(w)) or the Mellin transform Ai(s), it is known that the 
zeros/poles of these functions lie on certain lattices, and when (a,p) satisfy (7) these lattices overlap 
and almost all zeros/poles are cancelled (see Lemma 2 and the discussion after Corollary 2 in [25]). This 
explains why the general formula for the Wiener-Hopf factor (see (4.11) in [25]) or the general formula 
(22) for the Mellin transform can be simplified when X e Ck,i- The spectrally-positive { spectrally- 
negative } processes correspond to Doney's class Co,i { resp. C-1,-1 }, thus it seems that the classes 
Ck,i should be considered as generalizations of spectrally one-sided processes. An important problem is 
to a give probabilistic interpretation of Doney's classes. Hopefully, an answer to this question would 
also provide a probabilistic and more insightful derivation of the existing results on the Wiener-Hopf 
factorization, Mellin transform and the density of supremum in the case when X 6 Ck,i- 

Problem 2. In Theorem 2 we have established that the double series (12), (13) do not converge 
absolutely by showing that there is an infinite subsequence of the terms of the series which is bounded 
away from zero. However, this double series could still converge conditionally to p(x). Is it possible to 
find a way to order the partial sums of the series (12), (13) so that they converge to p(x) for all irrational 



Problem 3. From Theorems 1 and 2 we know that when a £ C U Q the series (12), (13) converge 
absolutely for all p, and when a G C these series do not converge absolutely for almost all p. This 
situation is reminiscent of Kolmogorov's zero-one law. Can this fact be established rigorously? Let us 
state this question more precisely: does there exist a value of a for which the series (12), (13) converge 
absolutely on a set of p of positive Lebesgue measure, and at the same time do not converge absolutely 
on a set of p of positive Lebesgue measure? 
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